package scu.maqiang.electromagnetics;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

import java.util.HashSet;
import java.util.stream.IntStream;

public class Scattering_2D {
    public static void main(String[] args) {
        //FEM2D_scat();
        FEM2D_scatDiel();
        //FEM2D_scat_layer();
        //FEM_scat_aniso();
        //FEM_scat_aniso2();
    }

    public static void FEM2D_scat() {

        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 300;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

        double mur = 1; // relative permittivity of object (used if scat_type='diel')
        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')


        /*
        //导入网格资源
        String meshFile = "RingMesh.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        //mesh.toTecplot("RingMesh.dat");

        //设置PML区域编号
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> {
            double r = Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r > 1.5 && r < 2.0)?1.0:-1.0;
        };
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        //设置与PEC接触边界
        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        //设置计算区域外边界
        int outerBoundaryLabel = 0;
        //ScalarFunc outerBounary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) > 1.5?1.0:-1.0;
        //mesh.setBoundaryLabel(outerBounary, null, outerBoundaryLabel);

        //设置散射区域与PML区域边界
        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 1.5) < 1.0e-3?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_Circle1PEC.msh");

         */

        /*
        String meshFile = "RectangleMesh_Circle.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0])< 1.5 && Math.abs(xy[1]) < 1.5?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;
//        ScalarFunc outerBounary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) > 1.5?1.0:-1.0;
//        mesh.setBoundaryLabel(outerBounary, null, outerBoundaryLabel);

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[1]) < 1.5?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[0]) < 1.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_Circle2PEC.msh");

         */


        /*
        String meshFile = "RectangleMesh_HalfCircle.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 1.5 && Math.abs(xy[1] - 0.5) < 1?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> {
            if ((Math.abs(xy[1]) < 1.0e-5 && Math.abs(xy[0]) < 1) ||
                    (Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.2 && xy[1] > 0.0) ){
                return 1.0;
            }
            return -1.0;
        };
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[1] - 0.5) < 1.0?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] - 0.5) - 1.0) < 1.0e-5) {
                return Math.abs(xy[0]) < 1.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_HalfCirclePEC.msh");
        */

        /*
        String meshFile = "RectangleMesh_Ellipse.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 1.5 && Math.abs(xy[1]) < 1?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.2?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[1]) < 1.0?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1]) - 1.0) < 1.0e-5) {
                return Math.abs(xy[0]) < 1.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_EllipsePEC.msh");

         */

        /*
        String meshFile = "RectangleMesh_LShape.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0] - 0.25) < 1.75 && Math.abs(xy[1] + 0.25) < 1.75?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> xy[0] > -1.5 && xy[0] < 2  && xy[1] > -2 && xy[1] < 1.5 ?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0] - 0.25) - 1.75) < 1.0e-5) {
                return Math.abs(xy[1] + 0.25) < 1.75?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] + 0.25) - 1.75) < 1.0e-5) {
                return Math.abs(xy[0] - 0.25) < 1.75?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_LShapePEC.msh");

         */


        /*
        String meshFile = "RectangleMesh_Conecircle.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 2.1 && Math.abs(xy[1] -0.575) < 2.675?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.abs(xy[0]) < 2.1 && Math.abs(xy[1] -0.575) < 2.675?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);


        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 2.1) < 1.0e-5) {
                return Math.abs(xy[1] - 0.575) < 2.675?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] - 0.575) - 2.675) < 1.0e-5) {
                return Math.abs(xy[0]) < 2.1?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_ConecirclePEC.msh");

         */

        /*
        String meshFile = "RectangleMesh_Ogive.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 2.5 && Math.abs(xy[1]) < 0.95?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> -pmlDomain.action(xy, label, param);
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 2.5) < 1.0e-5) {
                return Math.abs(xy[1]) < 0.95?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1]) - 0.95) < 1.0e-5) {
                return Math.abs(xy[0]) < 2.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);

        mesh.saveMesh("Scaterring_OgivePEC.msh");
         */


        //导入计算网格

        //1. 圆形区域PEC散射问题网格
        String meshFile = "C:\\Users\\maqia\\OneDrive\\meshFile\\Scattering_Circle1PEC.msh";
        //2. 圆形区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_Circle2PEC\\Scattering_Circle2PEC.msh";
        //3. 圆形区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_HalfCirclePEC\\Scattering_HalfCirclePEC.msh";
        //3. 椭圆区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_HalfCirclePEC\\Scattering_HalfCirclePEC.msh";
        //4. 椭圆区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_EllipsePEC\\Scattering_EllipsePEC.msh";
        //5. L型区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_LShapePEC\\Scattering_LShapePEC.msh";
        //6. 锥圆体PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_ConecirclePEC\\Scattering_ConecirclePEC.msh";
        //7. 卵形体PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_OgivePEC\\Scattering_OgivePEC.msh";

        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(MeshFormat.msh, meshFile);

        int PMLDomainLabel = 10;
        int innerBoundaryLabel = 1;
        int outerBoundaryLabel = 0;
        int interfaceBoundaryLabel = 3;

/*        HashSet<Integer> bSet = mesh.extractBoundaryNodes(3);
//        System.out.println(bSet.size());
//        for (int idx:bSet) {
//            System.out.println(idx + 1);
//        }

//        HashSet<Integer> bdSet = mesh.extractBoundaryNodes(3);
//        Object[] bdObj = bdSet.toArray();
//        int[] bdArray = new int[bdObj.length];
//        for(int i = 0; i < bdObj.length; i++) {
//            bdArray[i] = (int)(bdObj[i]);
//        }
//        Arrays.sort(bdArray);
//        System.out.println(bdArray.length);
//        System.out.println(Arrays.toString(bdArray));

//        Pair<Integer, double[]>[] PMLnodes = mesh.getNodeCoordinates(PMLDomainLabel);
//        System.out.println(PMLnodes.length);
//        for(int i = 0; i < PMLnodes.length; i++) {
//            System.out.println(PMLnodes[i].getFirst() + "\t" + Arrays.toString(PMLnodes[i].getSecond()));
//        }

//        Map<Integer, Complex[]> pmlCoordinates = mesh.setPMLCoordinates(k0, new int[]{interfaceBoundaryLabel}, new int[]{outerBoundaryLabel}, PMLDomainLabel);
//        for(var entry: pmlCoordinates.entrySet()) {
//            System.out.print(entry.getKey() + "\t");
//            Complex[] xyz = entry.getValue();
//            for(int i = 0; i < xyz.length; i++) {
//                System.out.print(xyz[i] + "\t");
//            }
//            System.out.println();
//        }

 */

        //入射角
        double phii = 0;            // degree, angle of incident field
        double varphi = phii * Math.PI/180;  // radian, angle of incident field
        //入射波函数的负值
        ScalarCFunc incidentWaveFunc = (xy, label, param) -> {
            Complex result = new Complex(0.0, 0.0);
            result.real =  -Math.cos(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            result.image = -Math.sin(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            return result;
        };

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SCMatrix A_PML = new SCMatrix(fs.GetNdof());

        //在自由空间组装整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, 0);

        //在PML区域组装PML矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                               A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);
        //整体矩阵
        SCMatrix K = A_PML.add(1.0, A);

        //右端项
        Complex[] RHS = Complex.zeros(fs.GetNdof());

        //整个外边界设为0值
        fs.applyBC_MBN(K, Direct.All, outerBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), outerBoundaryLabel);

        //PEC边界，结构区域内部边界
        fs.applyBC_MBN(K, Direct.All, innerBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{incidentWaveFunc}, innerBoundaryLabel);


        IterSCSolver solver = new IterSCSolver(K);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);

        //入射场负值
        Complex[] uint = mesh.functionInMesh(incidentWaveFunc, null);

        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Scattering_Circle2PEC.dat", result);
    }

    //
    public static void FEM2D_scatDiel() {
        /*
        String meshFile = "RingMesh2.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        mesh.toTecplot("RingMesh2.dat");

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> {
            double r = Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r > 1.5 && r < 2.0)?1.0:-1.0;
        };
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int DielDomainLabel = 20;
        ScalarFunc dielDomain = (xy, label, param) -> {
            double r = Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r < 1.0)?1.0:-1.0;
        };
        mesh.setDomainLabel(dielDomain, null, DielDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 1.0) < 1.0e-3?1.0:-1.0;
        mesh.addInterfaceBoundary(innerBoundary, null, innerBoundaryLabel);
        mesh.displayBoundarys();

        int outerBoundaryLabel = 2;
        ScalarFunc outerBounary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) > 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(outerBounary, null, outerBoundaryLabel);

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 1.5) < 1.0e-3?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        */


        /*
        //导入网格资源
        String meshFile = "RingMeshDiel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        //mesh.toTecplot("RingMesh.dat");

        String scattererElement = "ScattererEle_Circle1.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        //设置PML区域编号
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> {
            double r = Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r > 1.5 && r < 2.0)?1.0:-1.0;
        };
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        //设置与PEC接触边界
        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        //设置计算区域外边界
        int outerBoundaryLabel = 0;
        //ScalarFunc outerBounary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) > 1.5?1.0:-1.0;
        //mesh.setBoundaryLabel(outerBounary, null, outerBoundaryLabel);

        //设置散射区域与PML区域边界
        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 1.5) < 1.0e-3?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scattering_Circle1Diel.msh");

         */



        /*
        String meshFile = "RectangleMesh_Circle2Diel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);

        String scattererElement = "ScattererEle_Circle2.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0])< 1.5 && Math.abs(xy[1]) < 1.5?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;
//        ScalarFunc outerBounary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) > 1.5?1.0:-1.0;
//        mesh.setBoundaryLabel(outerBounary, null, outerBoundaryLabel);

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[1]) < 1.5?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[0]) < 1.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_Circle2Diel.msh");

         */

        /*
        String meshFile = "RectangleMesh_HalfcircleDiel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);

        String scattererElement = "ScattererEle_Halfcircle.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 1.5 && Math.abs(xy[1] - 0.5) < 1?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> {
            if ((Math.abs(xy[1]) < 1.0e-5 && Math.abs(xy[0]) < 1) ||
                    (Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.2 && xy[1] > 0.0) ){
                return 1.0;
            }
            return -1.0;
        };
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[1] - 0.5) < 1.0?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] - 0.5) - 1.0) < 1.0e-5) {
                return Math.abs(xy[0]) < 1.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scaterring_HalfcircleDiel.msh");

         */


        /*
        String meshFile = "RectangleMesh_EllipseDiel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);

        String scattererElement = "ScattererEle_Ellipse.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 1.5 && Math.abs(xy[1]) < 1?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.2?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 1.5) < 1.0e-5) {
                return Math.abs(xy[1]) < 1.0?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1]) - 1.0) < 1.0e-5) {
                return Math.abs(xy[0]) < 1.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scattering_EllipseDiel.msh");

         */

        /*
        String meshFile = "RectangleMesh_LShapeDiel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);
        String scattererElement = "ScattererEle_LShape.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0] - 0.25) < 1.75 && Math.abs(xy[1] + 0.25) < 1.75?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> xy[0] > -1.5 && xy[0] < 2  && xy[1] > -2 && xy[1] < 1.5 ?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0] - 0.25) - 1.75) < 1.0e-5) {
                return Math.abs(xy[1] + 0.25) < 1.75?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] + 0.25) - 1.75) < 1.0e-5) {
                return Math.abs(xy[0] - 0.25) < 1.75?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scattering_LShapeDiel.msh");
         */

        /*
        String meshFile = "RectangleMesh_ConecircleDiel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);

        String scattererElement = "ScattererEle_Conecircle.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 2.1 && Math.abs(xy[1] -0.575) < 2.675?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.abs(xy[0]) < 2.1 && Math.abs(xy[1] -0.575) < 2.675?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);


        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 2.1) < 1.0e-5) {
                return Math.abs(xy[1] - 0.575) < 2.675?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] - 0.575) - 2.675) < 1.0e-5) {
                return Math.abs(xy[0]) < 2.1?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scattering_ConecircleDiel.msh");
        */

        /*
        String meshFile = "RectangleMesh_OgiveDiel.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);

        String scattererElement = "ScattererEle_Ogive.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> Math.abs(xy[0]) < 2.5 && Math.abs(xy[1]) < 0.95?-1.0:1.0;
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> -pmlDomain.action(xy, label, param);
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        int outerBoundaryLabel = 0;

        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0]) - 2.5) < 1.0e-5) {
                return Math.abs(xy[1]) < 0.95?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1]) - 0.95) < 1.0e-5) {
                return Math.abs(xy[0]) < 2.5?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        mesh.saveMesh("Scattering_OgiveDiel.msh");

         */

        //mesh.displayBoundarys();
//        HashSet<Integer> bdSet = mesh.extractBoundaryNodes(3);
//        Object[] bdObj = bdSet.toArray();
//        int[] bdArray = new int[bdObj.length];
//        for(int i = 0; i < bdObj.length; i++) {
//            bdArray[i] = (int)(bdObj[i]);
//        }
//        Arrays.sort(bdArray);
//        System.out.println(bdArray.length);
//        System.out.println(Arrays.toString(bdArray));

//        Pair<Integer, double[]>[] PMLnodes = mesh.getNodeCoordinates(PMLDomainLabel);
//        System.out.println(PMLnodes.length);
//        for(int i = 0; i < PMLnodes.length; i++) {
//            System.out.println(PMLnodes[i].getFirst() + "\t" + Arrays.toString(PMLnodes[i].getSecond()));
//        }

        //1. 圆形区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_Circle1Diel\\Scattering_Circle1Diel.msh";
        //2. 圆形区域PEC散射问题网格
        //String meshFile = "D:\\MyJFem\\meshFile\\Scattering_Circle2Diel\\Scatterring_Circle2Diel.msh";
        //3. 圆形区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_HalfcircleDiel\\Scattering_HalfcircleDiel.msh";
        //4. 椭圆区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_EllipseDiel\\Scattering_EllipseDiel.msh";
        //5. L型区域PEC散射问题网格
        //String meshFile = "meshFile\\Scattering_LShapeDiel\\Scattering_LShapeDiel.msh";
        //6. 锥圆体PEC散射问题网格
        String meshFile = "meshFile\\Scattering_ConecircleDiel\\Scattering_ConecircleDiel.msh";
        //7. 卵形体PEC散射问题网格
//        String meshFile = "meshFile\\Scattering_OgiveDiel\\Scattering_OgiveDiel.msh";
//
//
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(MeshFormat.msh, meshFile);
        int DielDomainLabel = 20;
        int PMLDomainLabel = 10;
        int outerBoundaryLabel = 0;
        int interfaceBoundaryLabel = 3;


        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 300;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

//        Map<Integer, Complex[]> pmlCoordinates = mesh.setPMLCoordinates(k0, new int[]{interfaceBoundaryLabel}, new int[]{outerBoundaryLabel}, PMLDomainLabel);
//        for(var entry: pmlCoordinates.entrySet()) {
//            System.out.print(entry.getKey() + "\t");
//            Complex[] xyz = entry.getValue();
//            for(int i = 0; i < xyz.length; i++) {
//                System.out.print(xyz[i] + "\t");
//            }
//            System.out.println();
//        }

        //入射角
        double phii = 0;            // degree, angle of incident field
        double varphi = phii * Math.PI/180;  // radian, angle of incident field
        //入射波函数的负值
        ScalarCFunc incidentWaveFunc = (xy, label, param) -> {
            Complex result = new Complex(0.0, 0.0);
            result.real =  -Math.cos(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            result.image = -Math.sin(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            return result;
        };

        double mur = 1; // relative permittivity of object (used if scat_type='diel')
        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')


        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SCMatrix A_PML = new SCMatrix(fs.GetNdof());

        //组装介电区域左端矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, DielDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * 4}, BVPType.CONSIST_MASS, A, DielDomainLabel);

        Complex[] uint = mesh.functionInMesh(incidentWaveFunc, null);
        //组装右端源项1. 源项只在电介质内有效
        Complex[] RHS = Complex.zeros(fs.GetNdof());
        A.mul(uint, RHS);

//        for(int i = 0; i < fs.GetNdof(); i++) {
//            System.out.println(i + "\t" + RHS[i]);
//        }

        //在矩阵A上增加自由空间区域整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * 1.0}, BVPType.CONSIST_MASS, A, 0);

        //组装右端项2. 如果这样组装的话，最终会得到和书中一样的结果, 但觉得这样组装有问题, 因为在边界上会组装到在自由空间中的一部分内容
//        Complex[] RHS = Complex.zeros(fs.GetNdof());
//        A.mul(uint, RHS);
//        Complex[] newRHS = Complex.zeros(fs.GetNdof());
//        HashSet<Integer> nodeSet = mesh.extractDomainNodes(DielDomainLabel);
//        for(int idx:nodeSet) {
//            newRHS[idx] = RHS[idx];
//        }
//        RHS = newRHS;

        //组装PML区域整体矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);

        //得到所有区域整体矩阵
        SCMatrix K = A_PML.add(1.0, A);
        //K.sort();
        //System.out.println(K.display());
        //Complex[] RHS = Complex.zeros(fs.GetNdof());

//        ScalarCFunc innerBdFunc = (xy, label, param) -> {
//            Complex result = new Complex(0.0, 0.0);
//            result.real = -Math.cos(k0 * xy[0]);
//            result.image = -Math.sin(k0 * xy[0]);
//            return result;
//        };
//        Complex[] uint = mesh.functionInMesh(innerBdFunc, null);
//        Complex[] RHS = K.mul(uint);


//        fs.applyBC_MBN(K, Direct.All, outerBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), outerBoundaryLabel);
//
//        fs.applyBC_MBN(K, Direct.All, interfaceBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{innerBdFunc}, interfaceBoundaryLabel);

//        IterSCSolver solver = new IterSCSolver(K);
//        Complex[] x = Complex.zeros(fs.GetNdof());
//        solver.PCGSSSOR_2(RHS, x, 1.5, 2);
        DirectSSolver solver = new DirectSSolver(K);
        Complex[] x = solver.complexSolve(RHS, 1, 1e-6);

        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Scattering_ConecircleDiel.dat", result);
//        mesh.toTecplot("Scattering_2D2.dat", x);
    }


    public static void FEM2D_scat_layer() {
        double c0  = 3*1e8;           // m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double nmu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        //Input parameters
        //*****************
        double freq    = 3000;         // MHz, frequency
        freq    = freq*1e6;            // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0      = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg     = 2*Math.PI*freq;    // rad/sec, radial frequency
/*
        String meshFile = "RingMeshLayer.mmsh";
        Mesh2DT3Generator mg = new Mesh2DT3Generator();
        Mesh2DT3 mesh = mg.readMesh(meshFile, MeshFormat.mmsh);

        int layer1Label = 10;
        ScalarFunc layer1Region = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 0.11?1.0:-1.0;
        mesh.setDomainLabel(layer1Region, null, layer1Label);

        int layer2Label = 20;
        ScalarFunc layer2Region = (xy, label, param) -> {
          double r =  Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
          return (r > 0.11 && r < 0.12)?1.0:-1.0;
        };
        mesh.setDomainLabel(layer2Region, null, layer2Label);

        int PMLDomainLabel = 30;
        ScalarFunc pmlDomain = (xy, label, param) -> {
            double r =  Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r > 0.17)?1.0:-1.0;
        };
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);


        //设置与PEC接触边界
        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 0.15?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        //设置计算区域外边界
        int outerBoundaryLabel = 0;

        //设置散射区域与PML区域边界
        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 0.17) < 1.0e-5?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        //System.out.println(mesh.extractBoundaryNodes(interfaceBoundaryLabel).size());
        mesh.saveMesh("Scaterring_LayerPEC.msh");

 */

        String meshFile = "Scaterring_LayerPEC.msh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(MeshFormat.msh, meshFile);

        int layer1Label = 10;
        int layer2Label = 20;
        int PMLDomainLabel = 30;

        //设置与PEC接触边界
        int innerBoundaryLabel = 1;

        //设置计算区域外边界
        int outerBoundaryLabel = 0;

        //设置散射区域与PML区域边界
        int interfaceBoundaryLabel = 3;


        //入射角
        double phii = 0;            // degree, angle of incident field
        double varphi = phii * Math.PI/180;  // radian, angle of incident field
        //入射波函数的负值
        ScalarCFunc incidentWaveFunc = (xy, label, param) -> {
            Complex result = new Complex(0.0, 0.0);
            result.real =  -Math.cos(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            result.image = -Math.sin(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            return result;
        };

        double mur = 1; // relative permittivity of object (used if scat_type='diel')
        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')


        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SCMatrix A_PML = new SCMatrix(fs.GetNdof());

        //组装介电区域左端矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, layer1Label, layer2Label);
        ScalarFunc massCoeff = (xy, label, param) -> label == layer1Label?-k0 * k0 * 2:-k0 * k0 * 4;
        fs.assembleMass(new ScalarFunc[]{massCoeff}, null, BVPType.CONSIST_MASS, A, layer1Label, layer2Label);

        Complex[] uint = mesh.functionInMesh(incidentWaveFunc, null);
        //组装右端源项1. 源项只在电介质内有效
//        Complex[] RHS = Complex.zeros(fs.GetNdof());
//        A.mul(uint, RHS);

//        for(int i = 0; i < fs.GetNdof(); i++) {
//            System.out.println(i + "\t" + RHS[i]);
//        }

        //在矩阵A上增加自由空间区域整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * 1.0}, BVPType.CONSIST_MASS, A, 0);

        //组装右端项2. 如果这样组装的话，最终会得到和书中一样的结果, 但觉得这样组装有问题, 因为在边界上会组装到在自由空间中的一部分内容
        Complex[] RHS = Complex.zeros(fs.GetNdof());
        A.mul(uint, RHS);
        Complex[] newRHS = Complex.zeros(fs.GetNdof());
        HashSet<Integer> nodeSet = mesh.extractDomainNodes(layer1Label, layer2Label);
        for(int idx:nodeSet) {
            newRHS[idx] = RHS[idx];
        }
        RHS = newRHS;

        //组装PML区域整体矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);

        //得到所有区域整体矩阵
        SCMatrix K = A_PML.add(1.0, A);
        //K.sort();
        //System.out.println(K.display());
        //Complex[] RHS = Complex.zeros(fs.GetNdof());

//        ScalarCFunc innerBdFunc = (xy, label, param) -> {
//            Complex result = new Complex(0.0, 0.0);
//            result.real = -Math.cos(k0 * xy[0]);
//            result.image = -Math.sin(k0 * xy[0]);
//            return result;
//        };
//        Complex[] uint = mesh.functionInMesh(innerBdFunc, null);
//        Complex[] RHS = K.mul(uint);


        fs.applyBC_MBN(K, Direct.All, outerBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), outerBoundaryLabel);
        //整个外边界设为0值

        //PEC边界，结构区域内部边界
        fs.applyBC_MBN(K, Direct.All, innerBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{incidentWaveFunc}, innerBoundaryLabel);

        IterSCSolver solver = new IterSCSolver(K);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);

        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Scattering_Layer.dat", result);

    }

    public static void FEM_scat_aniso() {
        //1. 圆形区域PEC散射问题网格
//        String meshFile = "meshFile\\Scattering_Circle1Diel\\Scattering_Circle1Diel.msh";

        //导入网格资源
        String meshFile = "RingMeshAniso.mmsh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(meshFile, MeshFormat.mmsh);
        //mesh.toTecplot("RingMesh.dat");

        String scattererElement = "ScattererEle_CircleAniso.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        //设置PML区域编号
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> {
            double r = Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r > 1.5 && r < 2.0)?1.0:-1.0;
        };
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);

        //设置与PEC接触边界
        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        //设置计算区域外边界
        int outerBoundaryLabel = 0;

        //设置散射区域与PML区域边界
        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 1.5) < 1.0e-3?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(3).size());
        mesh.saveMesh(MeshFormat.msh, "Scattering_CircleAniso.msh");


        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 300;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

//        Map<Integer, Complex[]> pmlCoordinates = mesh.setPMLCoordinates(k0, new int[]{interfaceBoundaryLabel}, new int[]{outerBoundaryLabel}, PMLDomainLabel);
//        for(var entry: pmlCoordinates.entrySet()) {
//            System.out.print(entry.getKey() + "\t");
//            Complex[] xyz = entry.getValue();
//            for(int i = 0; i < xyz.length; i++) {
//                System.out.print(xyz[i] + "\t");
//            }
//            System.out.println();
//        }

        //入射角
        double phii = 0;            // degree, angle of incident field
        double varphi = phii * Math.PI/180;  // radian, angle of incident field
        //入射波函数的负值
        ScalarCFunc incidentWaveFunc = (xy, label, param) -> {
            Complex result = new Complex(0.0, 0.0);
            result.real =  -Math.cos(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            result.image = -Math.sin(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            return result;
        };

        double mur = 1; // relative permittivity of object (used if scat_type='diel')
        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')


        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SCMatrix A_PML = new SCMatrix(fs.GetNdof());

        //1. double[][] coef2 = {{1.0/9.0, 0}, {0, 1.0/4.0}};
        double[][] coef2 = {{1.0/4.0, 0}, {0, 1.0/9.0}};
        MatrixFunc mf = (xyz, label, param) -> coef2;
//        MatrixCFunc mf = (xy, label, param) -> {
//            Complex[][] coef = new Complex[2][2];
//            coef[0][0] = new Complex(4.0/15, 0);
//            coef[0][1] = new Complex(0.0, -1.0/15);
//            coef[1][0] = new Complex(0.0, 1.0/15);
//            coef[1][1] = new Complex(4.0/15, 0);
//            return coef;
//        };
        //装介电区域左端矩阵
        fs.assembleStiff(mf, null, BVPType.CONSTITUTE_MATRIX_COMMON, A, DielDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * mur}, BVPType.CONSIST_MASS, A, DielDomainLabel);

        Complex[] uint = mesh.functionInMesh(incidentWaveFunc, null);
        //组装右端源项1. 源项只在电介质内有效
//        Complex[] RHS = Complex.zeros(fs.GetNdof());
//        A.mul(uint, RHS);

//        for(int i = 0; i < fs.GetNdof(); i++) {
//            System.out.println(i + "\t" + RHS[i]);
//        }

        //在矩阵A上增加自由空间区域整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * mur}, BVPType.CONSIST_MASS, A, 0);

        //组装右端项2. 如果这样组装的话，最终会得到和书中一样的结果, 但觉得这样组装有问题, 因为在边界上会组装到在自由空间中的一部分内容
        Complex[] RHS = Complex.zeros(fs.GetNdof());
        A.mul(uint, RHS);
        Complex[] newRHS = Complex.zeros(fs.GetNdof());
        HashSet<Integer> nodeSet = mesh.extractDomainNodes(DielDomainLabel);
        for(int idx:nodeSet) {
            newRHS[idx] = RHS[idx];
        }
        RHS = newRHS;
        for(int i = 0; i < 5108; i++) {
            System.out.println(RHS[i]);
        }


        //组装PML区域整体矩阵
//        fs.assemblePMLStiffMass(mf, null, k0, -k0 * k0 * mur, BVPType.CONSTITUTE_MATRIX_COMMON, BVPType.CONSIST_MATRIX,
//                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);
        //得到所有区域整体矩阵
        SCMatrix K = A_PML.add(1.0, A);

//        K.sort();
//        System.out.println(K.display());
        //Complex[] RHS = Complex.zeros(fs.GetNdof());

//        ScalarCFunc innerBdFunc = (xy, label, param) -> {
//            Complex result = new Complex(0.0, 0.0);
//            result.real = -Math.cos(k0 * xy[0]);
//            result.image = -Math.sin(k0 * xy[0]);
//            return result;
//        };
//        Complex[] uint = mesh.functionInMesh(innerBdFunc, null);
//        Complex[] RHS = K.mul(uint);


//        fs.applyBC_MBN(K, Direct.All, outerBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), outerBoundaryLabel);
//
//        fs.applyBC_MBN(K, Direct.All, interfaceBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{innerBdFunc}, interfaceBoundaryLabel);

        IterSCSolver solver = new IterSCSolver(K);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);

        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Scattering_Aniso.dat", result);
        mesh.toTecplot("Scattering_2D2.dat", x);
    }


    public static void FEM_scat_aniso2() {
        //1. 圆形区域PEC散射问题网格
//        String meshFile = "meshFile\\Scattering_Circle1Diel\\Scattering_Circle1Diel.msh";

        //导入网格资源
        String meshFile = "RingMeshAniso.mmsh";
        Mesh2T3 mesh = new Mesh2T3();
        mesh.readMesh(meshFile, MeshFormat.mmsh);
        //mesh.toTecplot("RingMesh.dat");

        String scattererElement = "ScattererEle_CircleAniso.txt";
        int[] Scatter = MVO.readIntArray(scattererElement);
        int DielDomainLabel = 20;
        IntStream.of(Scatter).forEach(ele->mesh.setDomainLabel(ele - 1, DielDomainLabel));

        //设置PML区域编号
        int PMLDomainLabel = 10;
        ScalarFunc pmlDomain = (xy, label, param) -> {
            double r = Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]);
            return (r > 1.5 && r < 2.0)?1.0:-1.0;
        };
        mesh.setDomainLabel(pmlDomain, null, PMLDomainLabel);


        //设置与PEC接触边界
        int innerBoundaryLabel = 1;
        ScalarFunc innerBoundary = (xy, label, param) -> Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) < 1.5?1.0:-1.0;
        mesh.setBoundaryLabel(innerBoundary, null, innerBoundaryLabel);

        //设置计算区域外边界
        int outerBoundaryLabel = 0;

        //设置散射区域与PML区域边界
        int interfaceBoundaryLabel = 3;
        ScalarFunc interfaceBoundary = (xy, label, param) -> Math.abs(Math.sqrt(xy[0] * xy[0] + xy[1] * xy[1]) - 1.5) < 1.0e-3?1.0:-1.0;
        mesh.addInterfaceBoundary(interfaceBoundary, null, interfaceBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(3).size());
        mesh.saveMesh(MeshFormat.msh, "Scattering_CircleAniso.msh");


        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 300;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

//        Map<Integer, Complex[]> pmlCoordinates = mesh.setPMLCoordinates(k0, new int[]{interfaceBoundaryLabel}, new int[]{outerBoundaryLabel}, PMLDomainLabel);
//        for(var entry: pmlCoordinates.entrySet()) {
//            System.out.print(entry.getKey() + "\t");
//            Complex[] xyz = entry.getValue();
//            for(int i = 0; i < xyz.length; i++) {
//                System.out.print(xyz[i] + "\t");
//            }
//            System.out.println();
//        }

        //入射角
        double phii = 0;            // degree, angle of incident field
        double varphi = phii * Math.PI/180;  // radian, angle of incident field
        //入射波函数的负值
        ScalarCFunc incidentWaveFunc = (xy, label, param) -> {
            Complex result = new Complex(0.0, 0.0);
            result.real =  -Math.cos(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            result.image = -Math.sin(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            return result;
        };

        double mur = 1; // relative permittivity of object (used if scat_type='diel')
        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')


        FES2T31 fs = new FES2T31(mesh);
        SCMatrix A = new SCMatrix(fs.GetNdof());
        SCMatrix A_PML = new SCMatrix(fs.GetNdof());

        //1. double[][] coef2 = {{1.0/9.0, 0}, {0, 1.0/4.0}};
        //double[][] coef2 = {{1.0/4.0, 0}, {0, 1.0/9.0}};
        //MatrixFunc mf = (xyz, label, param) -> coef2;
        MatrixCFunc mf = (xy, label, param) -> {
            Complex[][] coef = new Complex[2][2];
            coef[0][0] = new Complex(4.0/15, 0);
            coef[0][1] = new Complex(0.0, 1.0/15);
            coef[1][0] = new Complex(0.0, -1.0/15);
            coef[1][1] = new Complex(4.0/15, 0);
            return coef;
        };
        //装介电区域左端矩阵
        fs.assembleStiff(mf, null, BVPType.CONSTITUTE_MATRIX_COMMON, A, DielDomainLabel);
        fs.assembleMass(new double[]{-k0 * k0 * mur}, BVPType.CONSIST_MASS, A, DielDomainLabel);

        Complex[] uint = mesh.functionInMesh(incidentWaveFunc, null);
        //组装右端源项1. 源项只在电介质内有效
//        Complex[] RHS = Complex.zeros(fs.GetNdof());
//        A.mul(uint, RHS);

//        for(int i = 0; i < fs.GetNdof(); i++) {
//            System.out.println(i + "\t" + RHS[i]);
//        }

        //在矩阵A上增加自由空间区域整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * mur}, BVPType.CONSIST_MASS, A, 0);

        //组装右端项2. 如果这样组装的话，最终会得到和书中一样的结果, 但觉得这样组装有问题, 因为在边界上会组装到在自由空间中的一部分内容
        Complex[] RHS = Complex.zeros(fs.GetNdof());
        A.mul(uint, RHS);
        Complex[] newRHS = Complex.zeros(fs.GetNdof());
        HashSet<Integer> nodeSet = mesh.extractDomainNodes(DielDomainLabel);
        for(int idx:nodeSet) {
            newRHS[idx] = RHS[idx];
        }
        RHS = newRHS;
        for(int i = 0; i < 5108; i++) {
            System.out.println(RHS[i]);
        }


        //组装PML区域整体矩阵
//        fs.assemblePMLStiffMass(mf, null, k0, -k0 * k0 * mur, BVPType.CONSTITUTE_MATRIX_COMMON, BVPType.CONSIST_MATRIX,
//                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);
        //得到所有区域整体矩阵
        SCMatrix K = A_PML.add(1.0, A);

//        K.sort();
//        System.out.println(K.display());
        //Complex[] RHS = Complex.zeros(fs.GetNdof());

//        ScalarCFunc innerBdFunc = (xy, label, param) -> {
//            Complex result = new Complex(0.0, 0.0);
//            result.real = -Math.cos(k0 * xy[0]);
//            result.image = -Math.sin(k0 * xy[0]);
//            return result;
//        };
//        Complex[] uint = mesh.functionInMesh(innerBdFunc, null);
//        Complex[] RHS = K.mul(uint);


//        fs.applyBC_MBN(K, Direct.All, outerBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), outerBoundaryLabel);
//
//        fs.applyBC_MBN(K, Direct.All, interfaceBoundaryLabel);
//        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{innerBdFunc}, interfaceBoundaryLabel);

        IterSCSolver solver = new IterSCSolver(K);
        Complex[] x = Complex.zeros(fs.GetNdof());
        solver.PCGSSSOR_2(RHS, x, 1.5, 2);

        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Scattering_Aniso2.dat", result);
    }

}
